# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)# charlibrary(stringr)
2 Setting up custom function
2.1windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",size = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",size =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",size = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",size =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")
4 Data Visualisation
4.1 Figure 1
Let’s create a table specific for this figure that only contains for each individual:
the first dive
the last dive
all benthic dives
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# then by individualgroup_by(id) %>%# keep only the first, last date, but also all benthic divesfilter(date ==min(date) | date ==max(date) | DiveTypeName =="Benthic")
Figure 1: Circular histogram plots displaying the times (in hours) of when female northern elephant seals (n = 403) perform their first recorded dive upon departing for their foraging trip (a), the last dive performed before returning from their foraging trip (b), and the time when they perform all benthic dives across their foraging trip (c).
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat)) %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName) %>%# create col to nicely display dive typemutate(origin =paste(DiveTypeName, "dives"))
Code
# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip <- trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1, 3, 2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# make it's group_by origindf_kernel_dens_sf <-group_by(df_kernel_dens_sf, origin)# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# facet by originfacet_wrap(. ~ origin)
Figure 2: Kernel density plots of the dives performed by female northern elephant seals (n = 403) during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Next step
see if it worth adding the North arrow and a scale
Using ggdensity instead of eks
This was the original way, but due to packages updates, I decided to switch with eks + ggOceanMaps option.
Code
# register to query google APIregister_google(key =readChar("../key/api_key.txt",file.info("../key/api_key.txt")$size))# get the maptrip_bis <-qmplot( long, lat,# weird, doesn't work if I increase nb of locations...data = df_kernel_dens %>%group_by(origin) %>%filter(row_number() <1000),geom ="blank",zoom =3,maptype ="satellite",source ="google")# mapstrip_bis +# kernelgeom_hdr(data = df_kernel_dens,aes(x = long, y = lat, fill =after_stat(probs)) ) +# remove some legendguides(alpha ="none") +labs(fill ="Probs") +# facetfacet_wrap2(. ~ origin)
Figure 3: Same but performed with ggdensity package
Figure 4: Overlap of female northern elephant seal’s trip (n=424) and predator habitat along with the percentage of benthic dives occurring inside as compared to outside of the habitat. The white star indicates the location of Año Nuevo State Park. Map a) displays orca and shark habitat along the west coast (adapted from Jorgenson, et al. 2019) and black lines represent the trip of all female northern elephant seals studied. Bar graph displays the average percentage of benthic dives occurring either inside or outside of orca habitat and displays the average percentage of benthic dives occurring inside and outside of shark habitat.
Figure 5: Overlap of female northern elephant seal’s trip (n=424) and predator habitat along with the percentage of benthic dives occurring inside as compared to outside of the habitat. The white star indicates the location of Año Nuevo State Park. Map a) displays orca and shark habitat along the west coast (adapted from Jorgenson, et al. 2019) and black lines represent the trip of all female northern elephant seals studied. Bar graph displays the average percentage of benthic dives occurring either inside or outside of orca habitat and displays the average percentage of benthic dives occurring inside and outside of shark habitat.
Next steps
see if it worth adding the North arrow and a scale
add the result of the test (inside vs. outside)
do we actually need legend for the map (area + bathymetry)?
4.4 Figure 4
Code
# main plotbathy_prop <- data_dive %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# create class of bathymetrymutate(bathy_class =fct_rev(cut( bathy,seq(-6000, 0, 400),ordered_result = T,dig.lab =4 ))) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per bathy_classgroup_by(bathy_class) %>%# to get the percentage of different dive types per bathy_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = bathy_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Bathymetry (m)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# second plotbathy_hist <- data_dive %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# create class of bathymetrymutate(bathy_class =fct_rev(cut( bathy,seq(-6000, 0, 400),ordered_result = T,dig.lab =4 ))) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = bathy_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotbath_plot <- bathy_hist / bathy_prop +plot_layout(heights =c(1, 10))
Code
# main plotdist_coast_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 1900, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the percentage of different dive types per dist_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from the coast (km)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# the second plotdist_coast_hist <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 2000, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = dist_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotdist_coast_plot <- dist_coast_hist / dist_coast_prop +plot_layout(heights =c(1, 10))
Code
# the main plotdist_dep_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_dep >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the percentage of different dive types per dist_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from departure (km)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# the second plotdist_dep_hist <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = dist_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotdist_dep_plot <- dist_dep_hist / dist_dep_prop +plot_layout(heights =c(1, 10))
Figure 6: Proportion of dives performed across entire trip separated by dive type across bathymetry intervals and distance from departure intervals, and displays the proportion of benthic dives occurring over the seals time spent at sea. For graph the largest percentage of benthic dives occurs at the lowest depth then occurs again at higher depths. For graph The highest percentage of benthic dives occurs at the lowest distance from where the seal is first assumed to have departed.
Next step
Try to keep a similar number of classes among figures
4.5 Figure 5
Code
# the last plotprop_at_sea <- data_dive %>%group_by(id) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# calculation of the proportionmutate(prop = nb_daily_dives_type / nb_daily_dives) %>%# order by seals/datearrange(id, nb_days_departure) %>%# perform our calculation per sealgroup_by(id) %>%# calculate of the percentage of time since departuremutate(perc_time_at_sea =round(nb_days_departure /max(nb_days_departure) *100, 1)) %>%# focus on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of percentage of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of percentage of time at seagroup_by(day_class) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="grey", col ="grey30") +guides(x =guide_axis(angle =45)) +scale_y_continuous(labels =function(x) { scales::percent(abs(x), 1) } ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Daily proportion of benthic dives (%)" )
Code
# displayprop_at_sea
Figure 7: Distribution of the percentage of benthic dives during elephant seals’ trip to sea expressed as a percentage of time spent at sea.
---title: "Figures"author: - name: "Astarte Brown" - name: "Joffrey Joumaa"date: "`r invisible(Sys.setlocale(locale = 'C')); format(Sys.Date(), format = '%B %d, %Y')`"format: html: toc: true toc-location: left number-sections: true smooth-scroll: true code-fold: true code-tools: true df-print: paged fig-align: "center" highlight-style: arrow self-contained: trueexecute: echo: true cache: false warning: falsetheme: light: flatly dark: darklyknitr: opts_chunk: message: false rownames.print: false tidy: styler---# Import library```{r}# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)# charlibrary(stringr)```# Setting up custom function## `windrose````{r}windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",size = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",size =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",size = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",size =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }```# Import dataLet's load `data_dive`, *i.e.* the output of `data_wrangling.qmd`, and filter only on animals leaving from Ano Nuevo.```{r}# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")```# Data Visualisation## Figure 1Let's create a table specific for this figure that only contains for each individual:* the first dive* the last dive* all benthic dives```{r}# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# then by individualgroup_by(id) %>%# keep only the first, last date, but also all benthic divesfilter(date ==min(date) | date ==max(date) | DiveTypeName =="Benthic")``````{r fig-windrose-1}#| fig-cap: "Circular histogram plots displaying the times (in hours) of when female northern elephant seals (n = 403) perform their first recorded dive upon departing for their foraging trip (a), the last dive performed before returning from their foraging trip (b), and the time when they perform all benthic dives across their foraging trip (c)."# for departurewindrose_departure <- data_windrose %>%group_by(id) %>%filter(DiveNumber ==1) %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Departure") %>%windrose(., grid =c(10, 20, 30), facet = T) +theme(axis.title.x =element_blank())# for arrivalwindrose_arrival <- data_windrose %>%group_by(id) %>%filter(date_tz ==max(date_tz)) %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Arrival") %>%windrose(., grid =c(10, 20, 30), legend_position ="top", facet = T) +theme(axis.title.y =element_blank())# for benthic diveswindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(3000, 6000, 9000), facet = T) +theme(axis.title.x =element_blank(),axis.title.y =element_blank() )# combine all(windrose_departure + windrose_arrival + windrose_benthic) +plot_layout(guides ="collect" ) &theme(legend.position ="top")```:::{.callout-caution}### Next step* Circular stats :):::## Figure 2:::{.callout-tip}This part is mostly based on [https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html](https://cran.r-project.org/web/packages/eks/vignettes/tidysf_kde.html). The [author](https://www.mvstat.net/tduong/) of the package confirms that his kernel density estimation is calculated using "meaningful" units such as UTM [^1], since the input is in a simple feature format.:::[^1]: * [https://gis.stackexchange.com/questions/64638/how-to-fill-in-the-parameters-for-kernel-density-estimation](https://gis.stackexchange.com/questions/64638/how-to-fill-in-the-parameters-for-kernel-density-estimation)* [https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1020&context=siv](https://digitalcommons.usf.edu/cgi/viewcontent.cgi?article=1020&context=siv)* [https://gis.stackexchange.com/questions/50692/choosing-projection-for-rasterization-of-random-latitude-and-longitude-data-in-n/50791#50791](https://gis.stackexchange.com/questions/50692/choosing-projection-for-rasterization-of-random-latitude-and-longitude-data-in-n/50791#50791)```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat)) %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName) %>%# create col to nicely display dive typemutate(origin =paste(DiveTypeName, "dives"))``````{r}# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip = trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1,3,2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}``````{r fig-map-1}#| fig-cap: "Kernel density plots of the dives performed by female northern elephant seals (n = 403) during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# make it's group_by origindf_kernel_dens_sf <-group_by(df_kernel_dens_sf, origin)# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# facet by originfacet_wrap(. ~ origin) ```:::{.callout-caution}### Next step* see if it worth adding the North arrow and a scale:::::: {.callout-tip collapse="true"}### Using `ggdensity` instead of `eks`This was the original way, but due to packages updates, I decided to switch with `eks` + `ggOceanMaps` option.```{r}#| eval: false#| echo: false# # to run if map need to be refreshed# # register token# register_google(key = "AIzaSyBRWtt6DPVlyAXQfQmxWO4QA4dSI6Vdoes")## # get the map from google# map_display_google_bw = get_googlemap(# center = c(lon = 180, lat = 30),# zoom = 2,# color = "bw")# map_display_google_terrain = get_googlemap(# center = c(lon = 180, lat = 30),# zoom = 2,# color = "color",# maptype = "terrain")# map_display_google_satellite = get_googlemap(# center = c(lon = 180, lat = 30),# zoom = 2,# color = "color",# maptype = "satellite")# map_display_google_satellite_bw = get_googlemap(# center = c(lon = 180, lat = 30),# zoom = 2,# color = "bw",# maptype = "satellite")# map_display_google_roadmap = get_googlemap(# center = c(lon = 180, lat = 30),# zoom = 2,# color = "color",# maptype = "roadmap")# map_display_stamen_terrain = get_map(# location = c(lon = -180, lat = 30),# zoom = 2,# color = "color",# source = "stamen",# maptype = "terrain")## # export# saveRDS(map_display_google_bw, "../export/map_display_google_bw.rds")# saveRDS(map_display_google_terrain, "../export/map_display_google_terrain.rds")# saveRDS(map_display_google_satellite, "../export/map_display_google_satellite.rds")# saveRDS(map_display_google_satellite_bw, "../export/map_display_google_satellite_bw.rds")# saveRDS(map_display_google_roadmap, "../export/map_display_google_roadmap.rds")# saveRDS(map_display_stamen_terrain, "../export/map_display_stamen_terrain.rds")# map_display_google_bw <- readRDS("../export/map_display_google_bw.rds")# map_display_google_terrain <- readRDS("../export/map_display_google_terrain.rds")# map_display_google_satellite <- readRDS("../export/map_display_google_satellite.rds")# map_display_google_satellite_bw <- readRDS("../export/map_display_google_satellite_bw.rds")# map_display_google_roadmap <- readRDS("../export/map_display_google_roadmap.rds")# map_display_stamen_terrain <- readRDS("../export/map_display_stamen_terrain.rds")``````{r fig-map-2}#| cache: true#| fig-cap: "Same but performed with ggdensity package"# register to query google APIregister_google(key =readChar("../key/api_key.txt",file.info("../key/api_key.txt")$size))# get the maptrip_bis <-qmplot( long, lat,# weird, doesn't work if I increase nb of locations...data = df_kernel_dens %>%group_by(origin) %>%filter(row_number() <1000),geom ="blank",zoom =3,maptype ="satellite",source ="google")# mapstrip_bis +# kernelgeom_hdr(data = df_kernel_dens,aes(x = long, y = lat, fill =after_stat(probs)) ) +# remove some legendguides(alpha ="none") +labs(fill ="Probs") +# facetfacet_wrap2(. ~ origin)```:::## Figure 3```{r}# get predator distribution (extraction using google earth)Shark_area_coord <-read_sf("../export/shark.kmz") %>%smooth(., method ="ksmooth", smoothness =0.5) %>%# st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%st_coordinates() %>%as_tibble() %>%select(lon = X, lat = Y)Orca_area_coord <-read_sf("../export/orca.kmz") %>%smooth(., method ="ksmooth", smoothness =0.5) %>%# st_crop(., xmin = 170, xmax = -110, ymin = 20, ymax = 59) %>%st_coordinates() %>%as_tibble() %>%select(lon = X, lat = Y)``````{r}# add inside/outside each areadata_dive <- data_dive %>%mutate(shark_area =point.in.polygon(point.x = Long,point.y = Lat,pol.x = Shark_area_coord$lon,pol.y = Shark_area_coord$lat ),orca_area =point.in.polygon(point.x = Long,point.y = Lat,pol.x = Orca_area_coord$lon,pol.y = Orca_area_coord$lat ),overlap_area =as.numeric(if_else(shark_area ==1& orca_area ==1, T, F)),full_area =as.numeric(if_else(shark_area ==0& orca_area ==0, F, T)) )``````{r}data_dive_longer <-pivot_longer( data_dive,c("shark_area","orca_area","overlap_area","full_area" ),names_to ="areas",values_to ="in_out")data_descriptive_ind <- data_dive_longer %>%filter(DiveTypeName =="Benthic") %>%group_by(areas, id) %>%summarise(n_tot =n(),n_inside =length(id[in_out ==1]),n_outside =length(id[in_out ==0]),.groups ="drop" ) %>%mutate(perc_inside = n_inside / n_tot,perc_outside = n_outside / n_tot )data_descriptive_area <- data_descriptive_ind %>%group_by(areas) %>%summarise(mean_perc_inside =wtd.mean(perc_inside, n_tot),mean_perc_outside =wtd.mean(perc_outside, n_tot),sd_perc_inside =sqrt(wtd.var(perc_inside, n_tot)),sd_perc_outside =sqrt(wtd.var(perc_outside, n_tot)),# based from here: https://stats.stackexchange.com/questions/25895/computing-standard-error-in-weighted-mean-estimationse_perc_inside =sqrt(wtd.var(perc_inside, n_tot)*sum((n_tot/sum(n_tot))^2)),se_perc_outside =sqrt(wtd.var(perc_outside, n_tot)*sum((n_tot/sum(n_tot))^2)), ) %>%pivot_longer(cols =ends_with("inside") |ends_with("outside"),names_to =c("stat", "in_out"),names_pattern ="(.*)_perc_(.*)" ) %>%pivot_wider(id_cols =c("areas", "in_out"),names_from ="stat",values_from ="value" ) %>%mutate(areas =factor(areas,level =c("overlap_area","orca_area","shark_area","full_area" ) ))data_stat_area <- data_dive_longer %>%filter(DiveTypeName =="Benthic") %>%group_by(areas, in_out) %>%summarise(nb_area =n(), .groups ="drop_last") %>%mutate(nb_total =sum(nb_area)) %>%summarise(rstatix::prop_test(x = nb_area, n = nb_total))``````{r}# geom_signif(# # y_position = data_descriptive_area %>%# # mutate(N = mean + sd) %>%# # select(N) %>%# # pull() %>%# # max() + 0.05,# y_position = 0.5,# xmin = c(0.6, 1.6, 2.6, 3.6),# xmax = c(1.4, 2.4, 3.4, 4.4),# annotations = if_else(# data_stat_area$p.signif == "****",# "***",# data_stat_area$p.signif# )# ) +``````{r}inside_bar_plot = data_descriptive_area %>%mutate(mean =if_else(in_out =="outside", -mean, mean),sd =if_else(in_out =="outside", -sd, sd),se =if_else(in_out =="outside", -se, se) ) %>%filter(in_out =="inside") %>%ggplot(aes(x = mean, y = areas, fill = areas)) +facet_grid(in_out ~ ., scales ="free_y") +geom_col(show.legend = F) +scale_fill_manual(values =c("#E1E826", "#E66100", "#9A6BEC", "#004D40")) +new_scale_fill() +geom_col(show.legend = F,mapping =aes(x = mean,y = areas,fill = in_out),alpha =0.3 ) +scale_fill_manual(values =c("grey70")) +geom_errorbar(mapping =aes(xmin =-se,xmax = se),# reduce size of horizontal bar of the errorbarwidth =0.05) +coord_flip(xlim =c(0, max(data_descriptive_area$mean))) +scale_x_continuous(breaks =c(-0.6,-0.4,-0.2, 0.2, 0.4, 0.6),expand =c(0, 0),labels =function(x) { scales::percent(abs(x), 1) } ) +scale_y_discrete(labels =function(x) {lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1]) } ) +theme(panel.spacing.y =unit(0, "mm"),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line.y =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),axis.title =element_blank(),axis.ticks.length.x =unit(0, units ="mm"),axis.text.x =element_blank() )outside_bar_plot = data_descriptive_area %>%mutate(mean =if_else(in_out =="outside",-mean, mean),sd =if_else(in_out =="outside",-sd, sd),se =if_else(in_out =="outside",-se, se) ) %>%filter(in_out =="outside") %>%ggplot(aes(x = mean, y = areas, fill = areas)) +facet_grid(in_out ~ ., scales ="free_y") +geom_col(show.legend = F) +scale_fill_manual(values =c("#E1E826", "#E66100", "#9A6BEC", "#004D40")) +new_scale_fill() +geom_col(show.legend = F, mapping =aes(x = mean, y = areas, fill = in_out),alpha =0.3) +scale_fill_manual(values =c("grey30")) +geom_errorbar(mapping =aes(xmin =-se,xmax = se ),# reduce size of horizontal bar of the errorbarwidth =0.05 ) +coord_flip(xlim =c(-max(data_descriptive_area$mean), 0)) +scale_x_continuous(# breaks = c(-0.6, -0.4, -0.2, 0.2, 0.4, 0.6),expand =c(0, 0),labels =function(x) { scales::percent(abs(x), 1) } ) +scale_y_discrete(labels =function(x) {lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1]) } ) +labs(y ="Areas") +theme(panel.spacing.y =unit(0, "mm"),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line.y =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed",ends ="first" )),axis.ticks.length.x =unit(0,units ="mm"),axis.title.y =element_blank())label_x_bar_plot =ggplot(data.frame(l ="Percentage of benthic dives", x =1, y =1)) +geom_text(aes(x, y, label = l, angle =90)) +theme(aspect.ratio =10,axis.title =element_blank(),axis.text =element_blank(),axis.ticks =element_blank()) +theme_void()``````{r}layout <-"ABBBBBBBBBBBACCCCCCCCCCC"predator_stat_choice_1 <- label_x_bar_plot +# (inside_bar_plot +theme(plot.margin =unit(c(0,0,0,0), "pt"))) + (outside_bar_plot +theme(plot.margin =unit(c(0,0,0,0), "pt"))) +plot_layout(design =layout)``````{r}# https://stackoverflow.com/questions/55015088/back-to-back-barplot-with-independent-axes-rpredator_stat_choice_2 <- data_descriptive_area %>%ggplot(aes(y = mean, x = areas, fill = areas)) +geom_bar(position ="fill", stat ="identity", show.legend = F) +scale_fill_manual(values =c("#E1E826", "#E66100", "#9A6BEC", "#004D40")) +new_scale_fill() +geom_bar(mapping =aes(y = mean,x = areas,fill = in_out),alpha =0.3,position ="fill", stat ="identity") +scale_fill_manual(values =c("grey70", "grey30")) +geom_errorbar(data = data_descriptive_area %>%filter(in_out =="outside"),mapping =aes(ymin = mean -se,ymax = mean + se ),# reduce size of horizontal bar of the errorbarwidth =0.05 ) +scale_y_continuous(labels =function(x) { scales::percent(abs(x), 1) } ) +scale_x_discrete(labels =function(x) {lapply(x, function(x) (str_split(str_to_title(x), "_"))[[1]][1]) } ) +labs(x ="Areas", y ="Percentage of benthic dives", fill ="Status") +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line.y =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),axis.ticks.length.x =unit(0,units ="mm"))``````{r}# color palette# https://davidmathlogic.com/colorblind/#%23D81B60-%231E88E5-%23FFC107-%23004D40# combine areasareas_coord <-rbind( Shark_area_coord %>%mutate(Predator ="shark"), Orca_area_coord %>%mutate(Predator ="orca"))# plotpredator_map <- trip +# new scalenew_scale_fill() +# add pointgeom_spatial_point(data = data_dive %>%filter(DiveTypeName =="Benthic"),aes(x = Long, y = Lat, fill ="black"),alpha =0.05,shape =20,stroke =NA,crs =4326 ) +# color palettescale_fill_manual(name ="Dive",labels ="Benthic",values ="black" ) +# should I display legend?guides(fill ="none") +# new scalenew_scale_fill() +# add spatial polygongeom_spatial_polygon(data = areas_coord,aes(x = lon,y = lat,fill = Predator,col = Predator ),alpha = .5,crs =4326 ) +# color palettescale_fill_manual(values =c("#E66100", "#9A6BEC") ) +scale_color_manual(values =c("#E66100", "#9A6BEC") ) +# set limit from predator_map$layers[[1]]$datacoord_sf(xlim=c(-5579139, 5579139), ylim =c(-8679599,-2642510))# reordering layers (continent on top)predator_map$layers <-c( predator_map$layers[c(1:2, 4:length(predator_map$layers))], predator_map$layers[3] )```:::{.panel-tabset}### Choice 1```{r fig-map-3}#| fig-cap: "Overlap of female northern elephant seal’s trip (n=424) and predator habitat along with the percentage of benthic dives occurring inside as compared to outside of the habitat. The white star indicates the location of Año Nuevo State Park. Map a) displays orca and shark habitat along the west coast (adapted from Jorgenson, et al. 2019) and black lines represent the trip of all female northern elephant seals studied. Bar graph displays the average percentage of benthic dives occurring either inside or outside of orca habitat and displays the average percentage of benthic dives occurring inside and outside of shark habitat."#| fig-height: 7# layoutlayout <-"#AAAAAAAAAAAAAAAAAAAAAAA#AAAAAAAAAAAAAAAAAAAAAAA#AAAAAAAAAAAAAAAAAAAAAAABCCCCCCCCCCCCCCCCCCCCCCCBDDDDDDDDDDDDDDDDDDDDDDD"# display(predator_map +theme(legend.position ="none")) + (label_x_bar_plot +theme(plot.margin =unit(c(0, 0, 0, 0), "pt"))) + (inside_bar_plot +theme(plot.margin =unit(c(0, 0, 0, 0), "pt"))) + (outside_bar_plot +theme(plot.margin =unit(c(0, 0, 0, 0), "pt"))) +plot_layout(design = layout)```### Choice 2```{r fig-map-3-bis}#| fig-cap: "Overlap of female northern elephant seal’s trip (n=424) and predator habitat along with the percentage of benthic dives occurring inside as compared to outside of the habitat. The white star indicates the location of Año Nuevo State Park. Map a) displays orca and shark habitat along the west coast (adapted from Jorgenson, et al. 2019) and black lines represent the trip of all female northern elephant seals studied. Bar graph displays the average percentage of benthic dives occurring either inside or outside of orca habitat and displays the average percentage of benthic dives occurring inside and outside of shark habitat."#| fig-height: 7# display(predator_map +theme(legend.position ="none")) / (predator_stat_choice_2 +theme(legend.position ="top")) +plot_layout(heights =c(1.5, 1))```::::::{.callout-caution}### Next steps* see if it worth adding the North arrow and a scale* add the result of the test (inside vs. outside)* do we actually need legend for the map (area + bathymetry)?:::## Figure 4```{r}# main plotbathy_prop <- data_dive %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# create class of bathymetrymutate(bathy_class =fct_rev(cut( bathy,seq(-6000, 0, 400),ordered_result = T,dig.lab =4 ))) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per bathy_classgroup_by(bathy_class) %>%# to get the percentage of different dive types per bathy_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = bathy_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Bathymetry (m)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# second plotbathy_hist <- data_dive %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# create class of bathymetrymutate(bathy_class =fct_rev(cut( bathy,seq(-6000, 0, 400),ordered_result = T,dig.lab =4 ))) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = bathy_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotbath_plot <- bathy_hist / bathy_prop +plot_layout(heights =c(1, 10))``````{r}# main plotdist_coast_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 1900, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the percentage of different dive types per dist_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from the coast (km)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# the second plotdist_coast_hist <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 2000, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = dist_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotdist_coast_plot <- dist_coast_hist / dist_coast_prop +plot_layout(heights =c(1, 10))``````{r}# the main plotdist_dep_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_dep >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the percentage of different dive types per dist_classmutate(percentage = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = percentage)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from departure (km)",y ="Dive type proportion (%)" ) +scale_fill_viridis("Dive Type",option ="plasma",discrete = T,direction =-1 ) +# scale_fill_manual("Dive Type", values = c("#fcfdbf", "#fc8961", "#b73779", "#51127c"))# position legend at the toptheme(legend.position ="none",panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) )# the second plotdist_dep_hist <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# ggplotggplot(aes(x = dist_class,y = N,fill = DiveTypeName,width =0.5 )) +geom_bar(stat ="identity",position ="dodge" ) +scale_fill_viridis(option ="plasma",discrete = T,direction =-1 ) +theme_void() +theme(legend.position ="top")# the actual plotdist_dep_plot <- dist_dep_hist / dist_dep_prop +plot_layout(heights =c(1, 10))``````{r fig-area-plot-4}#| fig-cap: "Proportion of dives performed across entire trip separated by dive type across bathymetry intervals and distance from departure intervals, and displays the proportion of benthic dives occurring over the seals time spent at sea. For graph the largest percentage of benthic dives occurs at the lowest depth then occurs again at higher depths. For graph The highest percentage of benthic dives occurs at the lowest distance from where the seal is first assumed to have departed."#| fig-height: 9# display(guide_area() + bath_plot + dist_coast_plot + dist_dep_plot) +plot_layout(nrow =5, guides ="collect")```:::{.callout-caution}### Next step* Try to keep a similar number of classes among figures:::## Figure 5```{r}# the last plotprop_at_sea <- data_dive %>%group_by(id) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# calculation of the proportionmutate(prop = nb_daily_dives_type / nb_daily_dives) %>%# order by seals/datearrange(id, nb_days_departure) %>%# perform our calculation per sealgroup_by(id) %>%# calculate of the percentage of time since departuremutate(perc_time_at_sea =round(nb_days_departure /max(nb_days_departure) *100, 1)) %>%# focus on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of percentage of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of percentage of time at seagroup_by(day_class) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="grey", col ="grey30") +guides(x =guide_axis(angle =45)) +scale_y_continuous(labels =function(x) { scales::percent(abs(x), 1) } ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Daily proportion of benthic dives (%)" )``````{r fig-prop-5}#| fig-cap: "Distribution of the percentage of benthic dives during elephant seals' trip to sea expressed as a percentage of time spent at sea."#| fig-height: 3# displayprop_at_sea```## Extra```{r}#| include: false#| eval: falsefwrite(data_dive[,.(id, date_start = date, date_end = date + Dduration, lat = Lat, lon = Long, dive_number = DiveNumber, benthic_index = BenthicDiveIndex, corner_index = CornerIndex, dive_type_name = DiveTypeName, dive_type = DiveType)], "./export/dive_data.csv")```